proc ds2;
  package math / overwrite = yes;
    method radians(double angle) returns double;
      return angle * constant('PI') / 180;
    end;

    method rational(double x, double maxden, in_out integer p, in_out integer q) returns double;
      /*
      ** FROM: https://www.ics.uci.edu/~eppstein/numth/frap.c
      ** FROM: https://stackoverflow.com/questions/95727/how-to-convert-floats-to-human-readable-fractions
      **
      ** find rational approximation to given real number
      ** David Eppstein / UC Irvine / 8 Aug 1993
      **
      ** With corrections from Arno Formella, May 2008
      **
      ** Modified for Proc DS2, Richard DeVenezia, Jan 2020.
      **
      ** usage: rational(r,d,p,q)
      **   x is real number to approx
      **   maxden is the maximum denominator allowed
      **   p is return for numerator
      **   q is return for denominator
      **   returns 0 if no problems
      **
      ** based on the theory of continued fractions
      ** if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...)))
      ** then best approximation is found by truncating this series
      ** (with some adjustments in the last term).
      **
      ** Note the fraction can be recovered as the first column of the matrix
      **  ( a1 1 ) ( a2 1 ) ( a3 1 ) ...
      **  ( 1  0 ) ( 1  0 ) ( 1  0 )
      ** Instead of keeping the sequence of continued fraction terms,
      ** we just keep the last partial product of these matrices.
      */

      declare integer m[0:1,0:1];
      declare double startx e1 e2;
      declare integer ai t result p1 q1 p2 q2;

      startx = x;

      /* initialize matrix */
      m[0,0] = 1; m[1,1] = 1;
      m[0,1] = 0; m[1,0] = 0;

      /* loop finding terms until denom gets too big */
      do while (1);

        ai = x;

        if not ( m[1,0] * ai + m[1,1] < maxden ) then leave;

        t = m[0,0] * ai + m[0,1];
        m[0,1] = m[0,0];
        m[0,0] = t;

        t = m[1,0] * ai + m[1,1];
        m[1,1] = m[1,0];
        m[1,0] = t;

        if x = ai then leave;     %* AF: division by zero;

        x = 1 / (x - ai);

        if x > 2147483647 /*x'7FFFFFFF'*/ then leave;  %* AF: representation failure;
      end;

      /* now remaining x is between 0 and 1/ai */
      /* approx as either 0 or 1/m where m is max that will fit in maxden */
      /* first try zero */

      p1 = m[0,0];
      q1 = m[1,0];
      e1 = startx - 1.0 * p1 / q1;

      /* now try other possibility */

      ai = (maxden - m[1,1]) / m[1,0];

      m[0,0] = m[0,0] * ai + m[0,1];
      m[1,0] = m[1,0] * ai + m[1,1];

      p2 = m[0,0];
      q2 = m[1,0];
      e2 = startx - 1.0 * p2 / q2;

      if abs(e1) <= abs(e2) then do;
        p = p1;
        q = q1;
      end;
      else do;
        p = p2;
        q = q2;
      end;

      return 0;
    end;
  endpackage;
run;
quit;


* Example;

proc ds2;
  data _null_;
    declare package math math();
    declare double x;
    declare int p1 q1 p q;

    method run();
      streaminit(12345);

      x = 0;
      do _n_ = 1 to 20;
        p1 = ceil(rand('uniform',9));
        q1 = ceil(rand('uniform',9));

        x + 1. * p1 / q1;

        math.rational (x, 10000, p, q);

        put 'add' p1 '/' q1 '  ' x=best16. 'is' p '/' q;
      end;
    end;
  enddata;
run;
quit;